function y=elasticity(beta,Scenario)
%This is a function to calculate elasticity using the estimates

N=size(beta,2);
S=size(beta,3);
L=size(beta,1);
option=L-size(Scenario,2)+1;

price_legal=Scenario(1,1);
price_street=Scenario(2,1);
dltime=Scenario(3,2);

beta_price=reshape(beta(4,:,:),N,S);
beta_dltime=reshape(beta(5,:,:),N,S);

w=[eye(option-1) Scenario];


prob=zeros(N,option-1);



for id=1:N
    elasl=zeros(1,3);
    elass=zeros(1,3);
    elasi=zeros(1,3);
    count=0;
    if S~=1
        for s=1:S
            count=count+1;
            numer=exp(w*beta(:,id,s));
            probi=(numer/(sum(numer)+1));

            for k=1:3
                elasl(k)=elasl(k)+beta_price(id,s)*price_legal*probi(1)*((k==1)-probi(k));
                elass(k)=elass(k)+beta_price(id,s)*price_street*probi(2)*((k==2)-probi(k));
                elasi(k)=elasi(k)+beta_dltime(id,s)*dltime*probi(3)*((k==3)-probi(k));
            end
            prob(id,:)=prob(id,:)+probi';
        end
    else
        for s=1:S
            count=count+1;
            numer=exp(w*beta(:,id,s));
            probi=(numer/(sum(numer)+1));

            for k=1:3
                elasl(k)=elasl(k)+beta_price(id,s)*price_legal*probi(1)*((k==1)-probi(k));
                elass(k)=elass(k)+beta_price(id,s)*price_street*probi(2)*((k==2)-probi(k));
                elasi(k)=elasi(k)+beta_dltime(id,s)*dltime*probi(3)*((k==3)-probi(k));
            end
            prob(id,:)=prob(id,:)+probi';
        end
    end
    prob(id,:)=prob(id,:)/count;
    elaslegal(id,:)=elasl/count;
    elasstreet(id,:)=elass/count;
    elasinternet(id,:)=elasi/count;
end

y=[mean(elaslegal./(ones(N,1)*mean(prob)));mean(elasstreet./(ones(N,1)*mean(prob)));mean(elasinternet./(ones(N,1)*mean(prob)))];